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Abstract 

In many time-dependent problems of practical interest the param- 
eters entering the equations describing the evolution of the various 
quantities exhibit uncertainty. One way to address the problem of 
how this uncertainty impacts the solution is to expand the solution 
using polynomial chaos expansions and obtain a system of differential 
equations for the evolution of the expansion coefficients. We present 
an application of the Mori-Zwanzig formalism to the problem of con- 
structing reduced models of such systems of differential equations. In 
particular, we construct reduced models for a subset of the polynomial 
chaos expansion coefficients that are needed for a full description of the 
uncertainty caused by the uncertain parameters. The viscous Burgers 
equation with uncertain viscosity parameter is used to illustrate the 
construction. For this example we provide a way to estimate the nec- 
essary parameters that appear in the reduced model without having to 
solve the full system. 



1 Introduction 

The problem of quantifying the uncertainty of the solution of systems of par- 
tial or ordinary differential equations has become in recent years a rather 
active area of research. The realization that more often than not, for prob- 
lems of practical interest, one is not able to determine the parameters, initial 
conditions, boundary conditions etc. to within high enough accuracy, has 
led to a flourishing literature of methods for quantifying the impact that this 
uncertainty imposes on the solution of the problems under investigation (see 
e.g. [61 HOI HH E21 HU US] ) . However, despite the increase in computational 
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power and the development of various techniques for uncertainty quantifica- 
tion there is still a wealth of problems where reliable uncertainty quantifica- 
tion is beyond reach. The main reason behind the inadequacy is the often 
high dimensionality (in probability space) of the uncertainty sources. When 
this uncertainty is coupled with the fact that for practical problems, even 
solving the corresponding equations for one value of the uncertain parameter 
(initial condition, boundary condition, . . .) can be very expensive, it results 
in the uncertainty quantification problem being a rather formidable task. 
One way to address this problem is to look for reduced models for a subset 
of the variables needed for a complete description of the uncertainty. 

We begin by noting that not all sources of uncertainty are created equal. 
For example, as we are taught from the theory of ordinary and partial differ- 
ential equations, the effect of uncertainty in the initial conditions is different 
from the effect of a parametric uncertainty (see e.g. [HE]). In addition, the 
effect of all types of uncertainty is intimately connected with the inherent 
instabilities that may be present in the underlying system which we sub- 
ject to the uncertainty. These considerations remain equally, if not more, 
important when we attempt to construct reduced models for uncertainty 
quantification. 

In the current work, we are concerned with the construction of reduced 
models for systems of differential equations that arise from polynomial chaos 
expansions of solutions of a PDE or ODE system. In particular, we focus 
on the case that the given PDE or ODE system contains some uncertain 
parameter and we want to construct a reduced model for the evolution of a 
subset of the polynomial chaos expansions that are needed for a complete 
description of the uncertainty caused by the uncertain parameters. There 
are different methods to construct reduced models for PDE or ODE systems 
(see e.g. 013] and references therein). We choose to use the Mori-Zwanzig 
(MZ) formalism in order to construct the reduced model [2 [3]. 

The main issue with all model reduction approaches is the computation 
of the memory caused by the process of eliminating variables from the given 
system (referred to as the full system from this point on) [3]. The memory 
terms are, in general, integral terms which account for the history of the 
variables that are not resolved. One would like, if possible, to compute these 
memory integrals without having to solve the full system. This is a difficult 
task, since it is rarely clear how the memory of a reduced model (which is 
based on the dynamics of the unresolved variables) can be estimated from 
pure analytical considerations or even relatively cheap numerical calculations 
involving only the resolved variables. On the other hand, for problems of 
practical interest where the solution of the full system may be, at best, only 
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feasible for short times, we are forced to consider ways of estimating the 
memory terms from such analytical or lower dimensional considerations. 

We use the case of the viscous Burgers equation with uncertain viscosity 
coefficient to illustrate how it is possible to estimate the parameters needed 
to specify the memory terms. The basic idea is that the uncertainty in 
the viscosity coefficient leads to linear and nonlinear contributions in the 
memory terms. One can group the linear contributions from the different 
memory terms and then require that this linear term is a stabilizing one. 
This procedure allows to estimate recursively (as we increase the order of the 
terms kept in the reduced model) the parameters involved in the memory 
integrals. 

Section [2] presents a brief introduction to the MZ formalism for the con- 
struction of reduced models of systems of ODEs. In Section [3] we develop a 
reformulation of the MZ formalism. This allows the calculation of the mem- 
ory terms through the solution of ordinary differential equations instead of 
the computation of convolution integrals as they appear in the original for- 
mulation. Section H] applies the reformulation of MZ presented in Section [3] 
to the viscous Burgers equation when the viscosity coefficient is uncertain. 
Finally, in Section [5] we discuss certain directions for future work. 



2 Mori-Zwanzig formalism 

We begin with a brief presentation of the Mori-Zwanzig formalism [2j [3]. 
Suppose we are given the system 

d ^-=R(t,u(t)), (l) 

where u = ({lifc}), k £ H U G with initial condition it(0) = Uq. Our goal is 
to construct a reduced model for the modes in the subset H. The system 
of ordinary differential equations we are given can be transformed into a 
system of linear partial differential equations 

^ = L0 fc , 4>k(uoi 0) = u ok , k e HUG (2) 

where L = Y^keHuG R ^ u o)d^-- Tne solution of © is given by u k (u ,t) = 
4>k (uo , t) . Using semigroup notation we can rewrite © as 

d 

-^e tL u ok = Le tL u ok 
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Suppose that the vector of initial conditions can be divided as no = (uo, no), 
where uq is the vector of the resolved variables (those in H) and Uq is the 
vector of the unresolved variables (those in G). Let P be an orthogonal 
projection on the space of functions of no and Q = I — P. 
Equation ^ can be rewritten as 

^-e tL u ok = e tL PLu ok + e tQL QLu ok + / e {t ~ s)L PLe sQL QLu ok ds, k £ H, 
dt J 

(3) 

where we have used Dyson's formula 



JL = e tQL + f e {ts)Lp Le sQL ds 

Jo 



(4) 



Equation ([3]) is the Mori-Zwanzig identity. Note that this relation is exact 
and is an alternative way of writing the original PDE. It is the starting 
point of our approximations. Of course, we have one such equation for each 
of the resolved variables u k ,k € H. The first term in ([3]) is usually called 
Markovian since it depends only on the values of the variables at the current 
instant, the second is called "noise" and the third "memory". 
If we write 

e tQL QLu ok = w k , 
Wk(uo,t) satisfies the equation 



§rw k (u ,t) = QLw k (u ,t) 



(5) 



[ w k (u , 0) = QLx k = R k (u ) - {PR k )(u ). 

If we project ([5]) we get 

d 

P—w k (u ,t) = PQLw k (u ,t) = 0, 

since PQ = 0. Also for the initial condition 

Pw k (u ,0) = PQLu ok = 

by the same argument. Thus, the solution of ([5]) is at all times orthogonal 
to the range of P. We call ([5]) the orthogonal dynamics equation. Since 
the solutions of the orthogonal dynamics equation remain orthogonal to the 
range of P, we can project the Mori-Zwanzig equation (J3j> and find 

^-Pe tL u ok = Pe tL PLu ok + P f e {t ~ s)L PLe sQL QLu ok ds. (6) 
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3 Finite memory 



In this section we describe a reformulation of the problem of computing the 
memory term which does not use the orthogonal dynamics equation. We fo- 
cus on the case when the memory has a finite extent only. The case of infinite 
memory is simpler and is a special case of the formulation presented below. 
Also, the current reformulation allows us to comment on what happens in 
the case when the memory is very short. 

Let w ok {t) = P f* e^ L PLe s ^ L QLu ok ds = P f* e sL PLe^- s ^ L QLu ok ds, 
by the change of variables t' = t — s. Note, that WQ k depends both on t and 
the resolved part of the initial conditions uq. We have suppressed the uq 
dependence for simplicity of notation. If the memory extends only for to 
units in the past (with to < t,) then 

Wok (t) = P [ e sL PLe it - s)QL QLu ok ds. 

J t-to 

The evolution of u>o k is given by 



dw ok 
dt 

where 



= Pe tL PLQLu ok - Pe {t - t(,)L PLe toQL QLu ok + w lk {t), (7) 



Wlk (t) = P f e sL PLe (t ~ s)QL QLQLu ok ds. 
J t-t 

To allow for more flexibility, let us assume that the integrand in the formula 
for Wi k (t) contributes only for t\ units with t\ < t . Then 

Wlk (t) = P f e sL PLe (t ~ s)QL QLQLu ok ds. 
Jt-ti 

We can proceed and write an equation for the evolution of w\ k (t) which 
reads 

= Pe tL PLQLQLu ok - Pe^"* 1 ) 1 PLe tlQL QLQLu ok + w 2k (t), (8) 



dt 
where 



w 2k (t) = P I e sL PLe^ QL QLQLQLu ok ds. 
Jt-u 



>t-ti 

Similarly, if this integral extends only for t 2 units in the past with t 2 < t\ , 
then 

W2k (t) = P [ e sL PLe ( - t '~ s)QL QLQLQLuo k ds. 

Jt-t-2. 
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This hierarchy of equations continues indefinitely. Also, we can assume for 
more flexibility that at every level of the hierarchy we allow the interval of 
integration for the integral term to extend to fewer or the same units of 
time than the integral in the previous level. If we keep, say, n terms in this 



Note that the last term in ([9]) involves the unknown evolution operator for 
the orthogonal dynamics equation. This situation is the well-known closure 
problem. We can stop the hierarchy at the nth term by assuming that 

Wnk{t) = 0. 

In addition to the closure problem, the unknown evolution operator for 
the orthogonal dynamics equation appears in the equations for the evolu- 
tion of w 0k (t), • • • , W(„-i)fc(i) through the terms P e ( *~* o)L PLe toQL QLu ok , . . . 
p e (t-to)Lp Le t QL(Q L -jn-iQ LuQk reS p ec tively. 

We describe now a way to express these terms involving the unknown 
orthogonal dynamics operator through known quantities so that we obtain 
a closed system for the evolution of lOofe(t), . . . , wr n _uf.(t). 

Since we want to treat the case where to is not necessarily small, we 
divide the interval [t — to,t] in no subintervals. Define 



where noAto = *o an d wok(t) = Ya=1 w ok Similarly, we can define the 




0) 



where 
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quantities w^(t), . . . , Wy^\t) 



.(!) 



(t) = P f e sL PLe {t - s)QL QLQLu ok ds 
Jt-Ati 



/•t-Afri 

P / e sL PLe {t ~ s)QL QLQLu ok ds 
Jt-2Ati 



w ^\t) = P e sL PLe (t - s)QL QLQLu ok ds, 



ft— (ni— l)Ati 
't-tl 

where niAti = ii and w\ k {t) = Y17=l w ik$)- I 11 a similar fashion we can 
define corresponding quantities for all the memory terms up to wr n _i\ k (t) = 

In order to proceed we need to make an approximation for the integrals 
over the subintervals. 

3.1 Trapezoidal rule approximation 

We have 



(i) 



W 0k 



(t)=P f e sL PLe {t - s)QL QLu Qk ds 
Jt-At 

= Pe tL PLQLu ok + Pe {t - Ato)L PLe AtoQL QLu ok 
from which we find 

Pe (t^t )L pLe At QL QLuok = ^_l_\ w (l)^ _ Pe tLp LQLuok + O((At ) 2 ) 

and from ([7]) 

dw$ ( 2 



~dt = ~ \, Al~ ) W ^ {t) + 2PetLpL Q Lu °x + <k (*) + O((A* ) a ). 
Similarly, for w^ k *(t) we find 

"Ut J ofei j 



2 



, , ,w$(t) - 2Pe tL PLQLuQ k + w ( $(t) + 0((At Q ^ 
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In general, 



dw. 



W 
Ok 



dt 

-i-l 



+ 



Ate 



4fc (*) + {-iy +1 2Pe tL PLQLu 0k 



E(^)(-ir +i+1 -S(*) 



+ W g(t) + 0((Ato) 2 ) for i = l,..., no- 



Similarly, 



(10) 




lk 



dt 

i-l 



Ah 



«>$(*) + (-l) l+ '2Pe tL PLQLQLu 0k 



+ 42W + 0((Ati) 2 ) for * = l,...,m 



(n-l)fc 



( l 1)fc W + (-l) l+1 2Pe tL PL(QL)™- 1 QL Uofe 



+ 



i-l 



■j=l 



+ 0((Ai n _!) 2 ) for i = 1, 



,n n _i. 
(11) 



By dropping the 0((Ato) 2 ), • • • , 0((Ai n _i) 2 ) terms we obtain a system of 
no + n\ + . . . + n n _i differential equations for the evolution of the quan- 



tities w, 



(i)/ 



(n-i)fc' This system allows us to determine the memory 
term wok(t) = P J e^~ s ^ L PLe s ® L QLuokds. Since the approximation we 
have used for the integral leads to an error O(At) 2 , the ODE solver should 
also be O(Ai) 2 . We have used the modified Euler method to solve numeri- 
cally the equations for the reduced model. 

Note that the implementation of the above scheme requires the knowl- 
edge of the expressions for Pe tL PLQLu k, ■■■ , P^ L PL{QL) a ^ x QLu^ k . Since 
the computation of these expressions for large n can be rather involved for 
nonlinear systems (see Section H]), we expect that the above scheme will 
be used with a small to moderate value of n. Finally, we mention that the 
above construction can be carried out for integration rules of higher order 
e.g. Simpson's rule. 



S 



4 Burgers equation with uncertain viscosity coef- 
ficient 



In this section we show how the above MZ formulation can be used for uncer- 
tainty quantification (UQ). In particular, we apply it to the one-dimensional 
Burgers equation when the viscosity coefficient is uncertain. The equation 
is given by 

u t + uu x = vu xx , (12) 

where v > 0. Equation (|12h should be supplemented with an initial condition 
u(x, 0) = uo(x) and boundary conditions. We solve (fT2|) in the interval [0, 2ir] 
with periodic boundary conditions. This allows us to expand the solution 
in Fourier series 

u N (x,t) = ^2u k (t)e ikx , 

where F = {— y , y — 1]. The equation of motion for the Fourier mode u k 
becomes 

du k ik s-^ 2 

— = -— 2^ u pUq - vk u k . (13) 

p+g=k 
p,qGF 

We assume that the viscosity coefficient v is uncertain (random) and can 
be expanded as v(£) = uq + a£ where £ is uniformly distributed in [—1,1]. 
In the numerical experiments we have taken = 0.1 and a = 0.07. This 
means that the viscosity coefficient is allowed to take values in the interval 
[.03, 1.07]. The choice of the range allows us to compute an accurate solution 
for any viscosity coefficient in the range without having to employ a large 
number of Fourier modes. 

To proceed we expand the solution Uk(t,£,) for k £ F in a polynomial 
chaos expansion using Legendre polynomials which are orthogonal in the 
interval [—1, 1]. In particular, we have that 

where Li(£) is the Legendre polynomial of order i. For each wavenumber k 
we expand the solution Uk(t,£,) of (]16p in Legendre polynomials and keep 
the first M polynomials 

M-l 

Mt,0 » Yl Mt)Li(0, where £ ~ U[-l, 1]. (14) 

i=0 
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Similarly, the viscosity coefficient can be written as v = Y2i=o u i^i(0 with 
v\ = a since Lo(£) = 1 and L\(£) = £. Substitution of (fH|) in (fT6|) . use of 
the expansion of the viscosity coefficient and of the orthogonality properties 
of the Legendre polynomials gives 

M-l M-l 

k2 U l U kmClmr (15) 

«=0 m=0 

for k G F and r = 0, . . . , M - 1. Also 

£?[L t (0L m (QL r (Q] 

where the expectation is taken with respect to the uniform density 
on [—1,1]. The expectation on the denominator of the expression for c\ mr 
is E[Lr(£)] = J*^ (£)^(f£ = 2^tl: wnne the expectation on the numer- 
ator can be computed accurately using Gaussian quadrature with Legen- 
dre nodes. The Legendre polynomial triple product integral defines a ten- 
sor which has the following sparsity pattern: E[Li(£)L m (£)L r (£)] = 0, if 
l + m<r or l + r<m or m + r<l or l + m + r = odd [8] . Due to this 
sparsity pattern, for a given value of M only about 1/4 of the M 3 tensor 
entries are different from zero. The sparsity pattern will be used below (see 
Section I4.2[) to facilitate the estimation of the length of the memory. 



, - ., M-l M-l 

dUkr{t) IK x - x - x - 

^ = ~~ ~2 2^1 2^/ U pl U qrnClmr 

1=0 m=0 p+q=k 

p.qeF 



4.1 MZ reduced model 

To conform with the Mori-Zwanzig formalism we set 

^ M-l M-l M-l M-l 

Rkr(u) = — ^ ] ^ 1 ^ ^ UplUqmClmr ~ k ^ ] ^ ] v lUk m Cl mr , 
1=0 m=0p+q=k 1=0 m=0 

p,q€F 

where u = {uk r } for k £ F and r = 0, . . . , M — 1. Thus, we have 

du kr 



dt 



Rkr(u) (16) 



for k £ F and r = 0, . . . , M — 1. We proceed by dividing the variables in 
resolved and unresolved. In particular, we consider as resolved the variables 
u = {u^} for k £ F and r = 0, . . . , A — 1, where A < M. Similarly, the 
unresolved variables are u = {uk r } for k £ F and r = A, ... ,M — 1. In 
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the notation of Section [2] we have H = F U (0, . . . , A — 1) and G = F U 
(A, . . . , M — 1). In other words, we resolve, for all the Fourier modes, only 
the first A of the Legendre expansion coefficients and we shall construct a 
reduced model for them. 

The system (TIE]) is supplemented by the initial condition uq = (uo,uq). 
We focus on initial conditions where the unresolved Fourier modes are set 
to zero, i.e. uq = (uo,0). We also define L by 



To construct a MZ reduced model we need to define a projection operator P. 
For a function h(uo) of all the variables, the projection operator we will use 
is defined by P(h(u)) = P(h(uo,uo)) = h(u ,0), i.e. it replaces the value 
of the unresolved variables uq in any function h(uo) by zero. Note that this 
choice of projection is consistent with the initial conditions we have chosen. 
Also, we define the Markovian term 



The Markovian term has the same functional form as the RHS of the full 
system but is restricted to a sum over only the first A Legendre expansion 
coefficients for each Fourier mode. 

4.2 Number of memory terms and memory length 
4.2.1 Number of memory terms 

We have to decide on the number of terms that will be used in the expansion 
of the memory as well as the length of the memory kept for each term (see 
Section [3]) . The fact that we are considering the case of uncertain viscosity 
coefficient becomes important in choosing how many terms to keep in the 
memory expansion and what the memory length should be for each term. 
To see this we need to compute the first few terms in the expansion. For 
the first two terms PLQLuokr and PLQLQLuokr we find 




M-l 
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PLQLuokr = 2 x 



M-l A-l 

r. 

l=A m=Op+q=k 
P,qeF 

M-l M-l 

mr 

1=0 m=A 



(17) 



and 



PLQLQLuQkr = 2 x 



-2 x 



"~2~ 



^ M-l A-l 

X] E X PLQLuopiUOgmCimr (18) 

M-l M-l 

— ^ X X PLu Opl PLu OqmClrr. 
l=A m=0 p+g=fe 

M-l M-l 

mClmr 

l=Q m=A 

For the sake of simplicity, we restrict attention to the case when A = 1, so 
that we resolve only the zeroth term in the Legendre expansion. The linear 
(viscous) term in (|15p contributes a linear destabilizing term in PLQLuoko 
and a linear stabilizing term in PLQLQLuoko- To show this, we use the fact 
that the Legendre expansion of the viscosity coefficient has only the zero and 
first components nonzero, the properties of the defined projection operator 
P and the sparsity of the Legendre polynomial triple product. Through 
straightforward but tedious algebra we find that the viscous term in ([T5]) 
contributes the term 



fc 4 ^CioiCii «OfcO 



in PLQLuoko an d the term 



-k ^icoiicioicno-uofco 



in PLQLQLuoko- Indeed, the contribution to PLQLu^o is destabilizing 
and the contribution to PLQLQLuoko is a stabilizing term. Note that these 
contributions correspond to terms of the form u xxxx and 

t^xxxxxx 

respectively 

in real space. 
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With more effort one can compute the contributions of the linear viscous 
term to the memory terms PLQLQLQLuQ k Q and PLQLQLQLQLuq^q. One 
finds that the viscous term contributes the term 

to PLQLQLQLuoko and the term 

-fc 10 [i/oi>iCiio3)iiCioi + Z-Wciiocoiicmcioicm 

+^O^fciioCi2lC 22ClOlCll2 + Z^lClloCl2lCll2CoilCi l]u fcO 

to PLQLQLQLQLuQko- 

The pattern of alternating destabilizing and stabilizing contributions of 
the viscous term to the memory terms continues for higher order terms. The 
conclusion from this pattern is that one needs to keep the memory terms in 
pairs in order to guarantee the stability of the reduced model. Also, since 
these contributions correspond to higher and higher spatial derivatives in 
real space we have to use only a few pairs otherwise the reduced model will 
become extremely stiff. In our numerical experiments we have used only the 
first pair of memory terms, namely PLQLuQko and PLQLQLuoko- 



4.2.2 Length of the memory 

We will use the information obtained in Section [4T27T] to estimate the length 
of the memory. Ideally, we would like to estimate the values of to an d t\ 
without having to solve the full system. That would make the construction 
of the reduced model efficient and applicable in cases where the solution of 
the full system is expensive (or possibly unknown). 

We focus again on the case when A = 1 and we assume that we use only 
one subinterval to discretize the time integrals, i.e. Aio = io an d Ati = t\. 
If we keep only the terms PLQLuQko and PLQLQLuoko for the memory, 
the reduced model reads 



duko _ , ;. 



dt 



e tL PLu oko + w oko {t) (19) 



= 2Pe tL PLQLu 0k0 - ^-w oko (t) + w lk0 (t) (20) 

= 2Pe tL PLQLQLu 0k0 - ^w lk0 (t) (21) 
at t\ 



13 



We can solve ([20]) and ([2T]) formally and substitute in (fl~9|) to get 



du k0 
dt 



e tL PLu oko + [ e~ Xo(t - s) 2Pe sL PLQLu oko ds 
Jo 

+ ( e - x °^ s) ( e~ Xl{s - r) 2Pe TL PLQLQLu kodTds, 
Jo Jo 



(22) 



where Ao = 2/to and Ai = 2/t±. 

The quantities u>ofco(£) and wifco(i) have on the RHS of their equations 
of evolution the terms Pe tL PLQLuq^q and Pe tL PLQLQLuoko respectively. 
As we have seen, the linear contributions of the viscous term to PLQLuQko 
and PLQLQLuQko correspond to higher spatial derivatives. Due to the 
presence of k 4 and — k G in the linear terms in the expressions for PLQLiiQko 
and PLQLQLuQko respectively, those linear terms are going to have large 
value for large wavenumbers. At the same time, these linear terms are linear 
in Uko(t) which is expected to evolve more slowly. Thus, the linear terms in 
the expressions for PLQLuQko an d PLQLQLuq^q are expected to have (at 
least for large wavenumbers) large values and evolve slowly. As a result, we 
expect the quantities u>ofco(^) an d vJiko(t) to evolve faster than Uko{t). With 
this in mind, we expect the memory lengths to and t\ to be shorter compared 
to the time scale of evolution of Uko(t). The crudest approximation that one 
can make for the integrals in ()22[) are 

t i j. 

e ~Mts) 2e sL pLQLuokods ^ 2Pe tL PLQLuoko = -^2Pe tL PLQLu 0k0 
Ao 2 



and 

e -A (t- s ) f e -M s -^2e TL PLQLQLu oko ^ t -^2Pe tL PLQLQLuoko- 
o Jo 4 

These approximations for the integrals allow us to group together all the 
linear terms on the RHS of the equation for Uko{t). Indeed, putting together 
the linear contributions from the Markovian term and the two integral terms 
(after the approximation) we get the linear term 



j 2 ..,42 *o*i , 6 2 

-k v cooo + t k z^cioicno — k ^lCoiicroiCno 



u*o(t). (23) 



The expression in brackets in (|23p can be used to determine what should 
the values of to and t\ be so that the reduced model is linearly stable for all 
wave numbers k £ F. 
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The calculation of to and t\ is done in two steps. If we ignore the fe 4 
and /c 6 terms, then the only contribution is from the Markovian term and 
the bracketed expression is a parabola in k with negative values. If we also 
include the k 4 term the parabola changes into a double- well curve which can 
become greater than zero for some wave numbers depending on the value of 
to. In fact, we can estimate the minimum value of to for this to happen by 
solving the equation 

— kmax^OCOOO + C^mai"! ClOlCllO = 

where k max is the maximum wavenumber present in the solution. In fact, 
t nm = (i / oCooo)/(^max zy i c ioiCiio)- For t™" all the wavenumbers are linearly 
stable except for the wavenumber k max which is only marginally stable. 

Now, suppose that we set to equal to t™ in and we also include the k 6 
term. Then, the bracketed expression in (I23p becomes a 6th order negative 
curve. The addition of the negative definite /c 6 term provides us with an 
advantage. It allows us to increase to to values larger than t nm as long as 
ii is appropriately chosen to make sure that the all the wavenumbers are 
linearly stable. Of course, one should not increase to too much because a 
correspondingly large value of t\ in conjunction with the A: 6 factor can render 
the reduced model very stiff. Thus, the final criterion which allows us to 
determine to an d t\ uniquely is that, based on the linear stability domain 
of the numerical method, we pick to an d t\ so that the required step size 
for the reduced model is not smaller than the step size for the original (full) 
system. 

4.3 Numerical results 

In this section we present numerical results for the reduced model of the 
viscous Burgers equation with uncertain viscosity coefficient given by v = 
Y^i=o v iLi{€) with vo = 0.1 and v\ = 0.07. The solution of the full system 
was computed with N = 96 Fourier modes (F = [—48,47]) and the first 
7 Legendre polynomials (M = 7). The first 7 Legendre polynomials were 
enough to obtain converged statistics for the full system. The full system 
was solved with the modified Euler method with At = 0.001. 

The reduced model uses N = 96 Fourier modes but only the first Leg- 
endre polynomial, so A = 1. The memory length parameters in the reduced 
model were chosen to be to = 0.2 and t\ = 0.01632 according to the scheme 
presented in Section T4.2.2i In particular, this choice of memory length guar- 
antees linear stability of the reduced model when it is solved with the modi- 
fied Euler method with a step size of At = 0.001. We discretize the memory 
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integral with 1 subinterval, i.e. no = n\ = 1, Aio = 0.2 and Aii = 0.01632 
according to the notation of Section [3l With this choice of parameters the 
running time of the reduced model is about half of that of the full system. 

Note that according to our analysis in Section [3] there is a discrepancy in 
the local truncation error estimates of the trapezoidal rule and the modified 
Euler scheme. For the trapezoidal rule the local truncation error estimate 
is O((0.2) 2 ) for the memory term wok and O((0.01632) 2 ) for the memory 
term w\k- On the other hand, the local truncation error estimate for the 
modified Euler method is O((0.001) 2 ). To make the discrepancy disappear 
we must use more subintervals at the cost of making the reduced model 
evolution more expensive. We tried that but the accuracy of the results of 
the reduced model did not change. 

Figure [T] shows the evolution of the mean energy of the solution 

^) = ^2vrKo(t)| 2 

fcfG-F 

as computed from the full system (with M = 7 Legendre polynomials), the 
MZ reduced model with A = 1 without memory (keeping only the Markovian 
term) and the MZ reduced model with A = 1 with memory. The reduced 
model performs equally well with or without memory. 

Figure [2] shows the evolution of the mean squared I2 norm of the gradient 
of the solution 

G(t) = Y,^k 2 \u k0 (t)\ 2 

as computed from the full system (with M = 7 Legendre polynomials), the 
MZ reduced model with A = 1 without memory (keeping only the Markovian 
term) and the MZ reduced model with A = 1 with memory. It is obvious from 
the figures that the inclusion of the memory term improves considerably the 
performance of the reduced model. 

By looking at Figure [21 we see that the reduced model with memory pre- 
dicts a smaller value for the peak of G{t). We know that the term PLQLu$ko 
contributes a linear destabilizing term to the reduced model (see (|23l ). So, 
the obvious question to ask is if one can improve the accuracy of the reduced 
model with memory by increasing to and correspondingly t\ at the cost of 
making the evolution of the reduced model more expensive. As explained in 
Section 14.2.21 the increase in the cost will come from the increased stiffness 
of the reduced model. We have tried increasing to and ii and the results did 
not become more accurate. The reason for this lack of improvement is a sign 
that if one wishes to improve the accuracy of the reduced model, one needs 
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Figure 1: Evolution of the mean of the energy of the solution using only the 
firf ; T 11 




Time 



Figure 2: Evolution of the mean of the squared I2 norm of the gradient of 
the solution calculated using only the first Legendre polynomial. 
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to include higher order terms. In particular, one will have to add at least 
the next pair of terms, namely PLQLQLQLuQko an d PLQLQLQLQLuq^q 
(see discussion at the end of Section 14.2.11 as to why the terms need to be 
added in pairs). 

5 Discussion and future work 

We have presented the application of the Mori-Zwanzig formalism to the 
construction of reduced models for systems of differential equations resulting 
from polynomial chaos expansions of solutions of differential equations with 
parametric uncertainty. In particular, we presented a way that the reduced 
model can be reformulated so that instead of integro-differential equations 
one has to solve differential equations. The problem that arises in any 
reduced model with memory is to compute the length of the memory. If 
possible, one wishes to obtain the length of the memory without having to 
solve the full system. For the case of the viscous Burgers equation with 
uncertain viscosity coefficient, we presented a way to actually compute the 
length of the memory without having to solve the full system. Note that 
this construction readily applies also to other equations that include viscous 
dissipation e.g. the Navier-Stokes equations. 

The viscous Burgers example highlights two important issues that arise 
when one wants to construct reduced models for parametric uncertainty 
quantification. 

The first issue is how to estimate the length of the memory integrals 
when we keep in the reduced model more than one coefficient in the Legendre 
expansion, so that A > 1. In this case, the bracketed expression in (j23|) will 
be replaced by a A x A matrix. In order for the reduced model to be linearly 
stable, we have to require that the matrix is negative definite (or at least 
semidefinite) . Since the elements of the matrix depend on the quantities 
to, ti, . . . , we can use the negative definite restriction to estimate to, ti, . . . . 

The second issue is also related to the length of the memory but ad- 
dresses a different aspect. We have seen for viscous Burgers that because 
the memory terms correspond to higher derivatives in physical space, the 
lengths of the integrals for the different memory terms (in our example to 
and t\) decrease as the order of the memory term increases. This allowed 
us to use a crude short-memory approximation of the memory integral (see 
Section [4. 2.2j) . On the other hand, there are cases when the length of the 
integrals for the different memory terms can increase as the order of the 
memory term increases. In such cases the short-memory approximation of 
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the memory integrals will not work. 

A simple example which illustrates this behavior is that of a single de- 
caying linear ode 

du 

— = —ku 
dt 

where k ~ U[0, 1]. If one applies the procedure outlined in Sections [3] and 
13.11 it is easy to see that the terms PLQLuoko, PLQLQLuq^q, . . . decrease in 
amplitude and thus the required memory integrals lengths to, t\, . . . increase 
as we go to higher order terms. A crude approximation is to assume that 
all the memory kernels in (122ft (and for higher order terms) have the same 
decaying characteristic times, that is to = t\ = • • •)• Preliminary numerical 
calculations show that as we increase the order of the memory terms kept 
in the reduced model we also have to increase the value of to i n order to 
increase the accuracy. A detailed analysis will be presented elsewhere. 

When the uncertainty is due not to a parameter in the equations but 
due to the initial conditions, the criterion presented in Section 14.2.21 for 
selecting the length of the memory will not work. For example, in the viscous 
Burgers equation, where the viscous term is diagonal in Fourier space, if 
the viscosity has no uncertainty, the projection operator makes the viscous 
term part of the Markovian term. As a result, the viscous term will end up 
contributing in the memory terms but the corresponding contribution is a 
term which is nonlinear in the resolved variables. Thus, even if one groups 
the memory contributions from the viscous term, there is no simple linear 
stability criterion, like the one invoked in the current work, to facilitate the 
estimation of the memory length. The construction of reduced models for 
the case of uncertain initial conditions will be presented in a forthcoming 
publication [13]. Such a construction can also be applied to the problem of 
constructing reduced models for systems forced by random noise [S]. 

Finally, we mention that one can construct models which effect reduction 
both for the variables needed to describe uncertainty and the number of 
variables needed to describe the system for one realization of the uncertainty 
sources. This two-level reduction is imperative in situations where solving 
even for one realization of the uncertainty sources is very expensive e.g. 
atmospheric flows, fluid structure interactions. 
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